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ABSTRACT 

A fully Sinc-Galerkin method in both space and time is presented for fourth-order time- 
dependent partial differential equations with fixed and cantilever boundary conditions. The 
sine discretizations for the second-order temporal problem and the fourth-order spatial prob- 
lems are presented. Alternate formulations for variable parameter fourth-order problems are 
given which prove to be especially useful when applying the forward techniques of this 
paper to parameter recovery problems. The discrete system which corresponds to the time- 
dependent partial differential equations of interest are then formulated. Computational issues 
are discussed and a robust and efficient algorithm for solving the resulting matrix system is 
outlined. Numerical results which highlight the method are given for problems with both 
analytic and singular solutions as well as fixed and cantilever boundary conditions. 


1 Thi8 research was supported by the National Aeronautics and Space Administration under NASA Con- 
tract No. NAS1-18605 while the author was in residence at the Institute for Computer Applications in 
Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 
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1 Introduction 


The Sinc-Galerkin method for partial differential equations (PDE’s) has previously been 
developed for the model elliptic problem in two and three dimensions [1, 2], the parabolic 
problem in one and two dimensions [3, 2], and the second-order hyperbolic problem in one 
dimension [4]. The present work extends the method to fourth-order time- dependent prob- 
lems with various common boundary conditions. This extension is important for the very 
practical reason that the numerical solution of problems in this class is necessary in ap- 
plications ranging from the control of large flexible space structures to the development of 
robotics designs [5, 6, 7]. 

For clarity of development, the method will be presented for the linear fourth-order time- 
dependent problems 

£«(*,() s^(*,i) + ^(s/(x)0(x.())=/(x,i), 0 < x < 1 (>0 

u(0, t) = u(l,<) = 0, t > 0 (1.1) 

£«M)=|(M) = o, <>o 

du 

u(x, 0) = (x , 0) = 0, 0 < x < 1 

C/Z 

and 

Cu(x,t ) = f(x,t ), 0<x<l i > 0 

u(o,() = a(i), ( £/ |jr)(M) = 7«. ( >° 0- 2 ) 

£(»SsO<‘- , >-* *>• ‘ >0 

u(x,0) = -«-(*) 0) = 0, 0<*<1. 

QZ 

These formulations are generalizations of the equations which arise when using the Euler- 
Bernoulli theory to model beams with flexural rigidity EI(x) and fixed and cantilevered 
ends, respectively. The general 7 (t) and £(<) in (1.2) allow for the inclusion of boundary 
controllers. For ease of presentation, the boundary conditions in (1.1) and (1.2) will respec- 
tively be referred to as fixed and cantilever conditions throughout the paper. It is noted that 
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the methods of this work are easily extended to problems with simple and free boundary 
conditions with further details given in [8]. 

The construction of an approximate solution to problems of the form (1.1) or (1.2) com- 
monly begins with a Galerkin discretization of the spatial variable with time-dependent 
coefficients. This yields a system of ordinary differential equations which is solved via dif- 
ferencing techinques. Due to stability constraints on the discrete evolution operator, low 
order methods with small time steps are often required to obtain accurate approximations. 
In contrast, the method of this work implements a Galerkin scheme in time as well as space. 
Because the basis functions are tensor products of sine functions composed with suitable 
conformal maps, the method has the inherent advantage that the study of error analysis and 
matrix structure begins at the level of an ordinary differential equation. 

The fully Sinc-Galerkin method in space and time has many other salient features due 
both to the properties of the basis functions and to the manner in which the problem is 
discretized. First, the judicious choice of a conformal map provides approximate solutions 
to (1.1) and (1.2) which are valid on the infinite time interval rather than only on a trun- 
cated time domain. Furthermore, the optimal exponential convergence rate is maintained 
even in the presence of boundary singularities. Finally, the discrete system requires no nu- 
merical integrations to fill either the coefficient matrix or the right-hand side vector. All 
three features prove to be advantageous when solving both forward and inverse fourth-order 
problems. A drawback to the method is that it produces a full system in contrast to the 
banded matrices which are associated with finite difference and finite element methods. In 
part, the exponential convergence rate offsets this disadvantage. 

The foundations of the Sinc-Galerkin method are described in Section 2. The fundamental 
quadrature rule is given, and the exponential convergence rate of this method is stated. A 
thorough review of sine function properties can be found in [9] and [10]. 

In the next section, the discretization of the second-order temporal and fourth-order spa- 
tial operators is outlined. Two Bchemes for discretizing EI(x) are presented. These two 
schemes are motivated on the one hand by the forward problem and on the other hand by 
the inverse problem involving the recovery of EI(x ), given sampled data. In the first, EI( x) 
is differentiated directly, whereas in the second (as motivated by the parameter recovery 
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problem) EI(x ) is replaced by a finite dimensional expansion EI mm before differentiation. 
Attention is focussed on preserving the method’s exponential convergence rate while dis- 
cretizing EI(x) and adapting to varying boundary conditions. 

In Section 4, the one-dimensional results from the previous section are combined to yield 
methods which very accurately approximate the solutions to the fourth-order time-dependent 
problems (1.1) and (1.2). Two equivalent formulations for the resulting discrete system are 
presented and a very robust and accurate solution algorithm is outlined. 

Numerical results are presented in the fifth section. Of the many examples tested, those 
discussed in this section best exhibit the features necessary for the practical implementation 
of the Sinc-Galerkin method. The first and second examples illustrate the method as applied 
to problems with fixed boundary conditons while the third and fourth examples have can- 
tilever boundary conditions. The first and last examples have analytic solutions, the second 
example has an algebraic singularity, and the third example contains a logarithmic singular- 
ity. The numerical results demonstrate that the exponential convergence rate is maintained 
in all four cases. 


2 Sine Function Properties 


For the Sinc-Galerkin method, the basis functions are derived from the Whittaker cardinal 
(sine) function 


.smc(x) = 


sin 


(xx) 


— OO < X < oo 


7TX 


and its translates 


S(k,h)(x) = sine ^ ^ , h > 0. 

For h* = 2, three adjacent members of this sine family (S(k t /i*)(x), k = —1,0,1) are shown 
in Figure 1. 
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Figure 1. Three Adjacent Members (S(k,h*)(x),k = - 1 , 0 , 1 , h* = 2) of the Translated Sine 
Family. 

To construct basis functions on the intervals ( 0 , 1 ) and ( 0 ,oo), respectively, consider the 
conformal maps 

«*) = '» (7 -~) (2 1) 

and 

T(tw) = /n(to). (2.2) 

The map <f> carries the eye- shaped region 

De = j* = x + iy : Jar# ^ - - 

onto the infinite strip 

D S = {( = C-M, ; 1,1 <<<<!}. 

Similarly, the map T carries the infinite wedge 

Dw = | w = t + ia : |ar0(to)| < d < — 

1 2 

onto the strip £)$ . These regions are depicted in Figure 2. 


I<«f} 


4 



Figure 2. The Domains Ds, D E , and D\y. 

The sine gridpoints z k € (0,1) in D E will be denoted x k since they are real. Similarly, the 
gridpoints w k £ (0,oo) in D w will be denoted t k . Both are inverse images of the equispaced 
grid in Ds\ that is 

X k = ^{kh) = 

and 

t k = T~ l (kh) = e kh . 

To simplify notation throughout the remainder of this section, the pairs <f>, D E and T ,D W 
are referred to generically as D. It is understood that the subsequent definition and 
theorems hold in either setting. Furthermore, the inverse of x i 8 denoted by ip. 

The important class of functions for sine interpolation and quadrature is denoted B(D ) 
and defined next. 
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Definition 2.1. Let B(D) be the class of functions F which are analytic in D, satisfy 

/, r . \ F (z)dz\ “» 0, i -*± oo 

where L = : |a| < d < ^ j, and on the boundary of D (denoted dD) satisfy 

N(F) = L \H*)M < 

Jbd 


The following theorem for functions in B(D ) is found in [11]. 

Theorem 2.1. Let T be (0,1) or (0,oo) when % — 4> or T, respectively. If F £ B(D) and 
z i ~ 'P(jh') = X -1 0'*0>j = il, ±2, . . ., then for h > 0 sufficiently small 


L 


F ^- h ,kM 


< K x e 


— 2ir d/h 


(2.3) 


Theorem 2.1 illustrates the exponential convergence rate which is a trademark of sine meth- 
ods. There is a common occasion when it is possible to evaluate the infinite series appearing 
in (2.3), namely when integrating against S(k, h) o x ■ In general, however, the series must 
be truncated. With additional hypotheses, it is proven in [9] and [12] that the truncation 
need not be at the expense of the exponential convergence. 


Theorem 2.2. Assume F £ B(D ) and that there exist positive constants K,a and (3 such 
that 


Fill 

X’(t) 


< K 


exp(-a|x(r)|), t £ ip((—oo,0)) 
exp(-/?|x(T)|), t £ V'QO, oo)). 


Then for h sufficiently small 


(2.4) 


l m dz~h ft 


n*) 


- M 


(2.5) 


Theorems 2.1 and 2.2 are used to establish a uniform error bound when building an 
approximate solution to an ordinary differential equation (ODE). It should be noted that 
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the nature of the class B(D ) guarantees that the exponential convergence rate holds for 
many differential equations with singular solutions; that is, problems where the solution has 
an unbounded derivative on the boundary. By applying the scheme to select second- and 
fourth-order ODE’s, one can derive the fundamental matrices comprising the discrete sine 
system for the fourth-order time-dependent problems of interest. 


3 Sinc-Galerkin Systems for ODE’s 


In this section, the sine discretizations will be catalogued for the second-order temporal 
problem and two different fourth-order spatial problems distinguished by their boundary 
conditions. Alternate formulations for the variable-parameter fourth-order problems will 
also be given which prove to be especially useful when applying the forward techniques 
outlined in this paper to parameter recovery problems. 

In order to construct the discrete Sinc-Galerkin system for either the temporal or spatial 
problems, the following identities are needed. Let 



[S(p, h) o x (z)} 



1, i = p 

0 , i^P, 


(3.1) 



— 5(p, h)o X (z) 


Z-Zi 


0, 

( i~P )’ 


i = p 

i^P, 


W ^ h* 


JL 

d X 2 


S{p,h)o X (z) 


3 

“ (-2)(-iy-p 

I ( i-py ’ 


i — p 


(3.2) 


(3.3) 


= A 3 


^5 S(p,h)ox(z) 


= < 


0, 

{ (i - p) 3 


l —P 


[6-7T 2 (i-p) 2 ], i^p, 


(3.4) 
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and 


= h * ^ s (p. a )°xW 


d* 


= < 


*=** 


T’ 


» = p 


*^ p > 


(3.5) 


(* - P) 4 

denote the evaluation at the gridpoint of the sinc-map compositions and their derivatives 
with respect to the map %, 


3.1. The Spatial Problem: Fixed Boundary Conditions 

In [13], a thorough analysis of the Sinc-Galerkin method is given for linear fourth-order 
ODE’s with fixed boundary conditions. For purposes of constructing the sine discretization 
for (1.1), it suffices to review that procedure for 


Lu(x) = (E7(®)u ,/ (x)) w = /(*), 0 < x < 1 

u(0) = u(l) = 0 ( 3 - 6 ) 

tx'(O) = w'(lj = 0. 


Note that the interval (0,1) is for convenience only; adapting the map ^ (see (2.1)) generalizes 
the method to any finite interval (a, b). 

To define the Sinc-Galerkin approximation to (3.6), select the basis where 

5i(x) = S(t, h s ) o <f>(x) and take the approximate solution to be 

JV. 

«m.(x) = 53 u i s i( x )> m x = M x + N x + 1. (3.7) 

i=-M, 

The unknown coefficients {ttj} in (3.7) are determined by orthogonalizing the residual 
Lu mu - f with respect to the functions This yields the discrete system 


(LUm. ~ f , S p ) = 0 

for p = —M x , ■ » • , N x . The weighted inner product (•, •) is taken to be 

(F, G) = f F(x)G(x)w(x) dx 
Jo 


where 


tu(z) = — 
v ; (At 


W*))» 


(3.8) 


(3.9) 
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For a further discussion concerning the choice of weight, see [13]. 

Before invoking the quadrature rules, integration by parts is used to transfer the differ- 
entiation of u onto S p w, thus yielding the system 

f u(x)[EI(x)(S p (x)w(x))"]"dx = f f(x)S p (x)w(x)dx (3.10) 

Jo Jo 

for p = — M x , ■ ■ • , N x . With the weight choice (3.9), the boundary terms 

{(EIu")'{S p w) - (EIu")(Spw)' + u'{EI(Spw)'y - u(EI(S p w)')"}(x) ' (3.11) 

o 

vanish for essentially all problems of interest. 

Two approaches are distinguished by the treatment of the first integral in (3.10). In the 
traditional scheme, EI(x) is differentiated directly and the resulting integrals are approx- 
imated via sine quadrature rules. This scheme is direct and suitable for a large class of 
forward problems. The second, alternative, approach is motivated by the parameter recov- 
ery problem and differs from the first in that EI{x) is replaced by a sine expansion, EI mm , 
before quadrature is applied. Both approaches then proceed in the same manner whereby 
the system is expanded and the resulting integrals are evaluated via Theorem 2.2, or when 
possible, Theorem 2.1. 

The careful choice of the decay parameters a and (3 in (2.4) provides a means of balancing 
the asymptotic errors resulting from the quadrature and hence minimizes the system size. 
With regard to (2.4), the condition 

{ x“ 35 € fo, — ^ 

" ( 3 - 12 ) 

(1 

guarantees the decay needed to truncate the infinite quadrature rule. A less general but 
more convenient assumption than (3.12) is 

\EI(x)u{x)\ < Kx a +$( 1 - xf f t. (3.13) 


With a and (3 specified and M x chosen, the parameter selections 



(3.14) 


9 


and 


N x = UM. + 1 


(3.15) 


balance the asymptotic quadrature errors to at least 0(e(~* daMm ^ ). This rate results from 
the presence of a sine function in the integral. In the above, [•] denotes the greatest integer 
function. Note that if is an integer, (3.15) can be replaced by the selection N x = 

The discrete system for (3.6), using the traditional approach, can then be formulated 
as follows. Let l^ l \l = 0, 1,2, 3, 4 denote the m, x m x matrices whose pi-th entry is o$ 
from (3.1) - (3.5) and let T>(rj) be the diagonal matrix with entries r}(x_M,), • • • , y{ x N m )- 

_ *4 

The vector of unknowns tZ = * ■ • , «w„] T is then related to the known vector / = 

[/(»_**.), ■•• > /(*J/.)] T by 

A X U — Z>((^') _| )/ 


where 

A x = 


i r /Wp(a 4 ) + ^/ (3) P(a 3 ) + ^/ (2) P(a 2 ) + ^-/ (1) P(a *) + / {0) P(a 0 ) 


The functions ai(x),t = 0, 1,2, 3, 4 are given by 


a 0 = 


(EIw) 

~ T 


"" n (El'w) w t ( EI"w )" 




in im inn 

a x = i(EIw) m + 6(EIwy'?j- + A{EJwy^rr + EIw?-j-* 

<p f " (p ' 

i/;; jl// 

- 6(£J'ti0'^7 - 2EI'w? J 7 + EI”(w )?- + 2(EI"w)' 

(b f (p (p* 


,m 7 


a 2 = 6(EIw)"0' + \2{EIw)'<f>" + AEIwf + 3(Elw) 

-6 (EI'w)$ - 6EIwf + EI"w<f>\ 

a 3 = A(EIw)'(<f>y + 6 EIw<f>’<f>" - 2 EI'w(<j>’) 3 




and 


a 4 = EIw(<f>y. 


(3.16) 


(3.17) 


(3.18) 


(3.19) 


(3.20) 

(3.21) 

(3.22) 


Further details concerning the derivation of the system (3.16) and a thorough spectral anal- 
ysis of the component matrices can be found in [13]. 
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As mentioned previously, the treatment of the first integral in (3.10) yields various pertur- 
bations of the method which are advantageous in certain applications. One such application 
is the parameter recovery problem where an integral part of most numerical schemes for 
solving that problem is an accurate forward solver. With this in mind, the alternative ap- 
proach mentioned above is implemented wherein the term EI{x ) in (3.10) is expanded as a 
linear combination of weighted sine functions with four Hermite-like algebraic terms. These 
terms are added to accommodate the potentially nonzero function and derivative values of 
El at x = 0 and x = 1. Specifically, this parameter basis is taken to be 

r 

+ l(®)> ^ = ~M X + 1 

M x ) = | t>f:(x)5*(i), -M x + 2<k<N x -2 ( 3 - 23 ) 

1>N.- i( x )> k = N x — 1 

Here S fc (x) = S(k,h x ) o and the basis weight v E is taken to be 

v E (x) = w(x) = [x(l — x)]$. (3.24) 

The algebraic boundary basis functions are given by 

= (1 - x ) 2 [2x + 1], 

&jv.-i( x ) = * 2 [2(1 — x) + 1], 

&_Af.( x ) = I ( 1 “ X f 

and 

b N.i x ) = -**(1 ~ ®)- 

The finite dimensional approximation of El then takes the form 

EI m .(x) = £ c k i> k {x). (3.25) 

k=-M, 
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The number of basis functions used in the expansion is chosen so as to guarantee a square 
coefficient matrix. This is done to simplify the implementation of the method when applied 
to the PDE (1.1) of interest. 

A quick note should be made concerning the choice of basis and the manner of expand- 
ing EI mm . The two derivative- interpolating boundary basis functions are added so that this 
expansion of EI mm is the same as that used with cantilever or free boundary conditions. The 
choice of (3.24) for basis weight is certainly sufficient and proves to be beneficial when in- 
corporating this forward scheme into a numerical method for solving the parameter recovery 
problem. 

The expansion (3.25) is substituted into (3.10) and the resulting integrals are evaluated 
via Theorem 2.2 or Theorem 2.1 when possible. The decay condition (2.4) equates to the 
condition 

„ f x a+ 3, xe(o, |) 

\£l{x)u(x)\ < K l K ' (3.26) 

[(l-xf+t, xe(±,l) 

where the “homogeneous” part of El is 

El{x) = S/(x)-E/(0)5_M. + i(x)- J E;/(l)fe J v.-i(x)-^/ , (0)^ - Mm (x)-Er(l)b Nm (x). (3.27) 

The arguments leading to this condition are analogous to those presented in th e seco n d-order 
case as described in [14]. Again, this may be replaced by the more stringent requirement 

|£Z(x)u(x)| < Kx a+ 1(1 - x^+f . 

As before, the asymptotic errors are balanced by choosing h x and N x as specified in (3.14) 
and (3.15). 

With d and / defined as before and El expanded, the system for (3.6) using this alter- 
native approach can be written as 

A x a=P((f)-?)/ (3.28) 

where 

A x = [*< 2 >P(p* (1 )) + 2$< 3 >I >(j^x>) + * (4) P(?v< o>)]. (3.29) 


12 


The notation V(p^t)),l = 0, 1,2 denotes the diagonal matrices containing the components 
of the vectors 

p*(0 = ¥Qc, 1 = 0,1,2 (3.30) 

where c = \c_M m , • • • , The matrices j = 2,3,4 and = 0,1,2 are defined 

componentwise by 

l$ (i) ]p< = ■^y( 5 P u; ) 0) ( x *') ( 3 - 31 ) 

and 

[*(0],, = (3.32) 

The notation on the right-hand sides of (3.31) and (3.32) indicates the j-th and £-th deriva- 
tives, respectively. 

To illustrate the dependence of = 2, 3, 4 and ^ l \£ = 0, 1,2 on previously defined 

matrices, the respective expansions are listed below. The diagonal matrices T> and the 
matrices I^\l = 0, 1,2, 3, 4 have sizes consistent with the following range of the indices i,p 
and k (—M x <i<N x , -M t <p<N t , -M« + 2 < k < N x - 2). From (3.31) it follows that 

$(*) = -L/(’)p(u^') + L^vpw') + I^V( w"), (3.33) 

h h x 


$(3) = l/( 3 )p (u;(^) 2 ) + ^j/ (J) ^( Zw'<f>' + 3w<j>") 

( 3 “" + + W T ) + ^ (t) 

$(4) = J_ I<*)v(w(<f>') 3 )+ ± /(3)p(4u; , (^) 3 4- 6 WT) 

^6tw V + 12 w'<f>" + 4t u<f>'" + 3tu^-^ 

+ t/' )v t 4 ”"' + + + w 9 ) 




(w™\ 

\T)' 


For 1 = 0, 1,2 the m, x m, matrices in (3.32) are given by 


*<-> = [a. i a. +1 ; i «U, i Si?.] 
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where 6^ = [^(x-m.), • • • , ^(xw.)] T for k = — M x , —M x + 1, N x — 1, and N x . Again, the 
superscript t indicates the £th derivative. The m x x (m x — 4) matrices B^ are 

1 9 (0) = V(v E )I {0) , (3.37) 

5 (1) = -^-V(v E <f>')I^ + V{y E )lW (3.38) 

h x 

and 

B (J) = pV (vstf) 2 ) I W - +-VM' + 2 »;*')/(■> + Z>K)/<“>. (3.39) 

The negative signs that appear in the definitions of B ^ and B ^ result from the transposing 
of jf 1 ). Again, it is noted that in (3.37) - (3.39), the m x x (m x — 4) matrices I^\l = 0, 1,2 
have components 6$ as defined in (3.1) - (3.3). 

Thus, the fourth-order spatial problem (3.6) can be solved in a variety of ways using 
the Sinc-Galerkin method. For standard forward problems, the system (3.16) is often the 
most convenient to formulate and solve. If the forward solver is part of a numerical routine 
for solving the parameter recovery problem, then (3.28) is more useful since El is replaced 
by its finite dimensional approximation. Both approaches yield solutions u m . which are 
exponentially convergent approximations to the solution u of (3.6). 

S.2. The Spatial Problem: Cantilever Boundary Conditions 

A second set of fourth-order boundary conditions arise when modeling beams that are 
fixed at one end and free at the other. To extend the Sinc-Galerkin method to problems 
with these cantilever boundary conditions, consider the ODE 

Lu(x) = (E I(x)u"(x))" = /(®), 0 < x < 1 

u(0) = a, (EIu")( 1) = 7 (3.40) 

u'(0) = ?, (EIu")'( 1) = 3. 


A Sinc-Galerkin method to approximate the solution of (3.40) can be developed as follows. 
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Define the set of basis functions by 


«») = 


B_ Mt _ 4 (x), 

i = —M x — 4 

B-m.~ 3 ( 1 ), 

II 

1 

£ 

1 

CO 

B-m.-i{^), 

II 

1 

1 

to 

B-m.-i(x), 

i = -Mg - 1 

v(x)Si(x), 

—M x <i < N x 

Bn.+i(x), 

i = N x + 1 

Bj/.+lfe), 

i = N x + 2 

Bn.+ s(®), 

i = N x + 3 

5 iV.+4(x), 

i = N x + 4. 


Here £,(®) = S(i, h x ) o <f>(x ) and the basis weight v(®) is taken to be 

v(x) = [®(1 — ®)] 3 . 

The boundary basis functions are 

B-m.- i(*) = (1 - ®) 4 [20® 3 + 10® a + 4® + 1], 

B Nm+l {x ) = ® 4 [20(1 - ®) 3 + 10(1 - ®) a + 4(1 - ®) + 1], 
B-m.- a(*) = ®(1 — ®) 4 [10® a + 4® + 1], 
5 jv. +2 (x) = -® 4 (1 - ®)[10(1 - ®) a + 4(1 - ®) + 1], 
B-v.-sM = x’(l - x) 4 [2* + i] , 

Bn.M x ) = X 4 (l - i) J [2(1 - x) + 1 


and 


b-m.-,(x) = ix 3 (i - x y 


B N.+t( x ) = ~^X 4 f 1 - x) 3 . 


(3.41) 


(3.42) 
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A brief note concerning the choice of basis is in order at this point. First, since 
h x ) o ^(®)] , l = 1, 2, ■ • ■ , are undefined at x = 0 and x = 1, some basis modifi- 
cations must be made when solving problems with nonzero boundary conditions (see also 
the definition of in (3.23)). It is tempting to use fewer algebraic boundary basis functions 
and the basis weight 

u(x) — ®(1 — ®) 3 , 

but in many problems this results in nonzero boundary terms when integrating by parts. By 
using the basis weight v(x) = [x(l — x)] 3 and a full complement of algebraic terms, this pitfall 
can be avoided. Furthermore, the basis {£»} as defined in (3.41) can be used for problems 
with free boundary conditions, thus providing consistency to the method. 

The approximate solution is then defined to be 

u m.( x ) = 1(®) + a(®) + 7Cw.+3(®) + 

+U_m.-3C-M.-3(z) + + u N.+l(N m +l(x) 

(3.43) 

ujv.+aC/v.+afa) 

N. 

+ 2 Ui£i(x) 

where 75,/?, 7, and 6 are known and the coefficients {ttj are unknown. The quantities 

^ s eT{ l) 7 

and 

1, EJ'(l) , 

- ei( 1) [£/(i)r 

are well-defined since EI(x) is assumed positive on [0, 1]. Note that with the definition (3.41) 
for the basis {&}> u *».( x ) satisfies the boundary conditions; that is, 


u m .(0) = <5, (£?/<.)( 1) = 7 

<,,(o) = 7 , (£/<.)'( i) = J. 

The m, + 4 unknown coefficients in (3.43) are determined by orthogonalizing the residual 
with respect to the set of sine functions {S p }”l± 2 Ma _ 2 . This Petrov- Galer kin approach is in 
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contrast to those Galerkin methods in which the residual is orthognalized with respect to 
the basis and is done to take advantage of the exponential accuracy of point evaluation in 
the quadratures. This yields the discrete system 

(Lu mm -f,S p ) = 0 (3.44) 

for p = —M t — 2, • • • , N x + 2. The inner product (•, •) is that defined in (3.8) with the weight 
in thiB case taken to be 

w(x) = 1. (3.45) 

The difference between the weight function in (3.45) and that given in (3.9) is due to the 
presence of the basis weight u(x) in the definition of the basis {£»}• If the definition 

u h( x ) = ]C “iCiC*) 

«=-M. 

is made, then (3.41) and (3.44) can be combined to yield 

{LB_M m ~3, Sp)u~M m ~3 + 4, S p )u_M m ~4 + (LB Na+i , S p )u Nm+ i 

+(LBf / m+ 2) S p )u Nm+3 + (Lu h ,S r ) (3.46) 

= ( 7 , S r ) 

for p = — M x — 2, • • • , N x + 2 and 

J(x) = /(x) - - JlLB_ Af._ 2 (x) - 7 LB Nm+3 (x) - SLB Nu+4 (x). 

Integration by parts is applied to ( Luh , S p ) thus yielding the integral 

/ Uh(x)[£/(x)(S' p (x)u;(x)) w ] ,, dx (3.47) 

Jo 

(compare to (3.10)). The weight choice (3.45) is sufficient for guaranteeing that the boundary 
term (3.11) vanishes with u replaced by u>,. As in Section 3.1 there are two approaches here. 
In the traditional approach EI(x) is differentiated directly and in the alternative approach 
El is simply replaced by the finite dimensional approximate EI mm . The remaining inner 
products in (3.46) are evaluated directly via Theorem 2.1. 
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Proceeding with the first approach, as indicated by (3.12), the choice of weight w directly 
affects the decay conditions dictated by (2.9) of Theorem 2.2. For the weight w{x ) = 1 the 
condition 

\EI(x)U(x)\ < x“ +3 (l - xf +3 (3.48) 

guarantees sufficient decay so that the asymptotic errors resulting from the quadrature can 
be balanced by choosing h x and N x as specified in (3.14) and (3.15). The term £/(x) denotes 
that part of the true solution which is approximated by uj, and is given by the formal change 
of variables 


U{x) — u(x) — o:B_Af._i(x) — /^ 0 _m.- 2 ( x ) ~ jB Nm+3 (x) — 6Bff m + ^{x) 
- U "(0)£L m .- 3 (x) - u'"(0)B_m._4(*) 

-u(l)B Na+ i(x) - u '( 1 )# jv .+ 2 ( x ). 


(3.49) 


The discrete system for (3.40) can then be formulated as follows. Let o_m.~ 3, d-M.-4, 
ajv.+i, #jv.+ 2 and / denote the (m x -f 4) x 1 vectors containing the product of and 
the approximations to the inner products (LB_m.- 3 , S p ), (LB_ Ma _ 4 , S p ), (LB_ Nm+ i, S p ), 
( LB Nm+7 ,S p ), and (J, S p ), respectively. Hence the p-th entries of the respective vectors are 

(Si) p = ( EIB i)"( x p)-J^j 


(i = -M x - 3, -M x - 4, N x + 1, N x + 2) and 


(Dp = 7(x P ) 


<f>'(x p y 


Furthermore, let A m denote the (m x + 4) x m x matrix which results from the expansion of 
the inner product ( Liih , S p ). For a*(x), ( = 0, 1 , 2, 3, 4 as defined in (3.18) - (3.22), the matrix 
A m is given by 


A m = + i/< 3 >I3(va,) + iiWPfva,) 

+ _L /(»!,(„„,) + ;(»)©( voo). 
h x 

Here I w ,i = 0, 1,2, 3, 4 are (m x + 4) x m x matrices whose pi-th entry is given by from 
(3.1) - (3.5), and v is defined in (3.42). 
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The discrete system for the determination of the unknown coefficients {«,•} is given by 

A x il=f 

where the (m, + 4) x (m x + 4) matrix A x is defined to be 

A x = [d_jvf > _4 : 3_Af._ 3 : A m : Hn.+i : a^. +2 ]. (3.50) 

Here tf is defined to be the (m x + 4) x 1 vector 

3= «-M.» • ' ' » u A^.> u W,+I> u W.+2] T (3.51) 

containing the unknowns. 

It is noted that the matrices A x as defined in (3.17) and (3.50) differ only in the presence 
of v in the diagonal multipliers and the addition of border vectors. Hence the method is easily 
adapted when the boundary conditions are changed. Moreover, the exponential convergence 
rate is maintained, thus preserving the accuracy of the method. 

With parameter recover in mind, it is again worthwhile to UBe the alternative approach 
to develop the discrete system which arises when EI{x) is replaced by the finite dimensional 
term EI mm as defined in (3.25). To simplify notation in the discussion which follows, let 

AT. 

EIh(x) = c k v B (x)S k (x) 

k=-M. 

and 

EI e (x ) = 1(*) + c_Af._a&_M.-2(z) + cjv.+i&Tv.+^a:) + Ctf. +3 &Ar. + j(a:) 

so that EI mm = Elh + EI C . Note that EI k and EI C simply denote the sine and algebraic 
components of the approximate parameter. 

Consider now the inner products found in the system (3.46). The expansion of (Lu k , S p ) 
proceeds exactly as before with EI(x ) simply being replaced by EI ma (x) in the integral 
(3.47). In the boundary inner products, ( LB { , S p ),i = — M x — 4, • ■ • , —M x — 1 and i = N x + 1, 
• • • , N x -f 4, expansion and integration by parts yields 

(LBi, S p ) = {£ EI h (x)B'/{x){S p w)"{x)dx + H T J 

+ f\EI c B[')"(x)S p {x)w{x)dx. 

Jo 
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The weight w(x ) = 1 (see (3.45)) is sufficient for guaranteeing that 

B t = {(EI h B?)'(S p w) - EI h B?(S p w)'}(x)\l = 0. 

The resulting integrals are evaluated via Theorem 2.2 or, when possible, Theorem 2.1. 

For the weight u/(x) = 1 , the decay condition is 

\Sl{x)U(x)\ < x a+3 (l - xf f3 

where U and El are defined in (3.49) and (3.27) (compare to (3.48)). Again, the asymptotic 
errors are balanced by choosing h x and N x as specified in (3.14) and (3.15). 

The matrix system corresponding to (3.40) may be formulated as follows. Let 'id*', 
and p^(<) be defined as they were in (3.31), (3.32), and (3.30), respectively (with ^ = 0,1,2 
and j = 2,3,4). Note that in the definitions now, the index ranges are — M x <i< N x , 

—M x — 2 < p < N x + 2, and —M x — 2 < k < N x + 2, and the (m x + 4) x 1 coefficient vector 
is now B = • j C iv.+ 2 ] T - Hence and p^i) have the sizes (m x + 4) x m x , 

m x x (m x + 4) and m x x 1, respectively. 

Furthermore, let 4" denote the (m x + 4) X m s matrix which is defined componentwise by 

and let c = [c_Af., • • • , cjv.] t . Finally, for i = —M x — 4, • ■ ■ , — M x — 1 and i = AF X +1, • ■ • , N x + 4, 
let cti denote the (m x + 4) x 1 vectors 

3i= 4"V{v e B”)c + V ^(£/ c £")"j r] • ! 

Here 1 is simply the m x x 1 vector of ones. 

With the unknown vector u given by (3.51), the discrete system can be written as I 

A x xL=f \ 

where i 

I 

A x = [2_m.-4 : B_ Na _ 3 : A m : ajy. +1 : a^. +2 ] (3.52) ! 

| 

and 

—f \ 1 - 7 ~ 7^+3 - SaN+4- 
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The (m, + 4) x m x submatrix A m is given by 

A m = [$< a >£>(t>£*< s >) + 2*V>V(vfax)) + * (4) t>(titf»«»)]. 

It should be noted that the coefficient matrix A x in (3.52) differs from that arising in 
the fixed boundary problem, (3.29), only in the presence of v in the diagonal multipliers and 
the addition of border vectors. This makes the method easily adaptable when changing the 
boundary conditions. Furthermore, the matrices and can be expanded in terms 

of fundamental matrices in a manner similar to that in (3.33) - (3.36), thus simplifying the 
implementation of the method. Finally, the exponential convergence rate of the method is 
maintained, thus preserving the method’s accuracy. 

With the techniques from this section, the implementation of the Sinc-Galerkin method 
for problems with free and simple boundary conditions can be accomplished in a manner that 
is completely analogous to that used for cantilever boundary conditions. Further details and 
examples of the Sinc-Galerkin method for problems with free and simple boundary conditions 
can be found in [8]. 

S.S. The Temporal Problem 

The last ODE to be considered is the initial value problem 


Pu(t) = ti(<) = f{i ) , 0 < t < oo 
■u(O) = u(0) = 0. 


(3.53) 


A Sinc-Galerkin method to approximate the solution of (3.53) can be developed in a manner 
similar to that of the preceding boundary value problems. Define the set of basis functions 

by 

S;(t) = S(j,h,)oT(t) 

where T : Dw —* Ds is given in (2.2). The approximate solution ti m( (i) is then defined by 


N t 


= X) u 3 S i( t )> rn t = Mt + N t + 1. 
i=-M, 


(3.54) 


The m t unknown coefficients {uj} in (3.54) are determined by orthogonalizing the residual 
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with respect to the set of sine functions which leads to the analysis of 


(/, s;) = (Pu,s') 

= (*,•?,*) 

for q = —M t , • ■ ■ , Nf The weighted inner product for (3.55) is defined to be 

(F, G) = f" F{t)G{t)w'(t)it, 

J 0 

and the weight is taken to be 


(3.55) 


»•(*) = ~n 


1 


y/W) 

for reasons that are discussed in [4]. As before, integration by parts is used to transfer the 
differentiation of u onto S*w*. To guarantee that boundary terms vanish, it is assumed that 

lim -j-ttc- = 0 . 

*-*o ln(t) 

The resulting integrals are then evaluated via the quadrature rules of Section 2 . With respect 
to (2.4), the condition 

n — . I <€(0,1) 

WOv^WI < * 

[ <-*, < E [1, oo) 

guarantees the boundedness necessary to truncated the infinite quadrature rule. With 7 and 
6 specified and chosen, the parameter selections 


h t = 


j 7T d 

jM t 


and 


N t = 


i M > + ‘] 


(3.56) 


balance the asymptotic errors in (2.5) to at least ). 

In many time-dependent PDE’s, it is reasonable to assume that the solution decays 
exponentially at infinity; that is, the solution satisfies 


K')v / t«l < * 


r, is (o,i) 

e~ St , t e [ 1 , 00 ) 
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or, more stringently, 


KOI < KP+U- H . (3.57) 

With this supposition, Lund [12] shows that the condition (3,56) can be replaced by 

N t = [L/ n (Jma) + l] . (3.58) 

The selection N t in (3.58) significantly reduces the size of the discrete system with no loss 
of accuracy. 

The discrete system for (3.53) can then be formulated as follows. Let I^\l = 0, 2 denote 
the nij x mt matrices whose qj-th entry is 6$ from (3.1) and (3.3), and let t] ) again be 
the diagonal matrix with entries • • • »^(tAr«). With the usual definitions for i 1 and / 

and the identity 

(t(0)-*{(t(())-»}" = - j. 

the system for the determination of the unknown coefficients {uj} is given by 

A t a = V{{T)-l)f (3.59) 


where 


A t = 


1 lW _ 1 / 
4 


w 


P((f)i). 


(3.60) 


Further details concerning the derivation of the system (3.59) can be found in [4]. 

Note that nonzero initial conditions can be handled in a manner analogous to that used 
for nonzero boundary conditions in the previous discussion. Rational initial basis functions 
are used to incorporate the initial behavior and this known contribution is then taken to the 
right-hand side of the resulting discrete system (see [8]). All other analysis is identical to 
that for the problem (3.53). 


4 Time-Dependent Problems 

This section details the Sinc-Galerkin method applied to fourth-order time-dependent PDE’s 
with fixed and cantilever boundary conditions. Since the choice of basis, test functions, and 
inner product are all straightforward extensions of those used to solve the ODE’s in Section 3, 
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the error analysis and system formulation follow directly from previously discussed results. 
Once a discrete system has been formulated, various options exist for solving the associated 
matrix equation. Two such algorithms are outlined and their relative merits for various 
problems are discussed. 

4.1. The Time-Dependent Problem: Fixed Boundary Conditions 
Consider the time-dependent problem 

3 7 ( 3 7 u \ 

0 - *) + fat l *)) = ^ x ’ *)> 0 < * < 1 t > 0 ... 

(4.1) 

Sij(x, t ) = Si(x)Sj(t) = S(i } h x ) o <f>(x)S(j , ht) O T(<), 

the approximate solution is defined by way of the tensor product expansion 

N. N t 

12 12 m x = M x + N x + 1, m t = M t + N t + 1. 

i=-Afi j——Mt 

The m x ■ mt unknown coefficients {ujj} are determined by orthogonalizing the residual with 
respect to the set of sine functions {S p (x)S*(t)}^ZM^"--'!N m - This yields the discrete Galerkin 
system 

(Cu mmmt -f,s p s;) = 0 

for p = —M x , • • ■ , N x and q = — M t , • • • , N t . The inner product (•, •) is taken to be 

(F,G)= f f F(x,t)G(x,t)w(x,t)dxdt (4-2) 

Jo Jo 

where 

w(x,t) — w(x)w*(t) = (^'(a:)) - $(T(f)) - i. 

The quadrature rules and one-dimensional results from Sections 3.1 and 3.3 can be used to 
determine the resulting matrix system. 


u(0, t) = u(l,t) = 0, <>0 

= t>0 

3u 

«(*,°) = -^-(^, 0 ) = 0 , 0 < X < 1 . 

Given the basis {S^} where 


24 



As before, equating asymptotic errors is fundamental to minimizing system size. When 
the decay conditions (3.13) or (3.26), and (3.57) are combined to yield 


\EI{x)u{x,t)\ < AV* + *(1 -xf+U^e- St 

(43) 

or 

\SI{x)u(x,t)\ < Kx a ^{ 1 - xf+U^e-'*, 

(4.4) 

then the choices 

i 

II 

, =* 
, 

(4.5) 

V aM x 

II 

>- 

H 

(4.6) 

N. = + lj , 

(4.7) 

Mt = \-M, + l! , 

(4.8) 

IL7 J1 

and 

Wi =[-1„(1m A ) + i] 

(4.9) 

for the stepsizes and summation limits balance the asymptotic errors. 

If one takes d = 


then the above choices yield an asymptotic error rate of order 0(e V *^) 


Given M x , N x , M t , N t and h = h a = h t as defined above, the discrete system for (4.1) is 

A x UCj + C x UAj = G. (4.10) 


Here 


and 



The m t x m t matrix A t is given by 


A t 


-L M -ij 

hi 4 


V ((t)») 


(4.11) 

(4.12) 


(4.13) 
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as shown in (3.60). Furthermore, m x x m t matrices U and F are defined componentwise by 

[U]ij = Uij 

and 

= )• 

It should be noted that the ordering of the coefficients tty in U mimics that used in most 
standard time-differencing schemes. This is a matter of convenience since the Sinc-Galerkin 
method is not bound by any specific ordering of the grid. 

The structure of the tti x x m x matrix A x depends on the scheme that is used to discretize 
EI(x). If the parameter is fully differentiated, then A x is given by (3.17). If, on the other 
hand, EI(x ) is approximated by a linear combination of sine and algebraic basis functions, 
then A x is given by (3.29). 

Various methods exist for solving the equation (4.10). Referred to as a generalized 
Sylvester equation (4.10) is algebraically equivalent (page 414 of [15]) to the system 

Au = {C t <g> A x -f A t ® C x } co(U) = co(G) (4-14) 

where the tensor or Kronecker product of an m x m matrix E with a p x q matrix H is 
defined by 

E® H = [eijH] mpXnq . 

—4 

The vector U = co(U ) is the concatenation of the m x x m t matrix U obtained by successively 
"stacking” the columns of t/, one upon another, to obtain an m x m< x 1 vector. 

The system (4.14) can be solved directly via any of the decomposition methods that are 
available for linear systems. Although this system is easily formulated, the fact that A is 
very large (m x • m t x m x ■ m t ) and not banded causes this method to be impractical in some 
problems. For more general fourth-order operators however, this may be the only method of 
choice. 

A second algorithm for solving (4.10) depends on the generalized Schur decomposition 
(page 396 of [16]). As guaranteed by the results of Moler and Stewart [17], there exist unitary 


26 



matrices Qi,Z i} Q 2 , and Z 2 such that 


Q\A X Z X = P 


Q\C x Zi = R 


Q 2 C t Z 2 = S 
Q 2 AfZ 2 = T 


where P,R,S, and T are upper triangular. If Y = Z[UZ 2 and G = Q\GQ 2 , then (4.10) 
transforms to 


PYT* + RYS* = C. 


By comparing the fc-th columns, one finds that 


n n 

p Y + R Y s ^yj = c * 

j=k i~k 


which yields 

n n 

{hkP + SkkR)Vk — c k — P XI ^kiUj — R Y (^-^) 

j=fc+i i=fc+l 

(for convenience, it is assumed that all matrices are n x n and indexed from 1 to n). With 
the assumption that the matrix ( ikkP + s kkR) * s nonsingular, the solution to (4.15) is easily 

found by recursively solving triangular systems. 

Although this algorithm does require complex algebra, it is both robust and efficient 
and requires no assumptions concerning the diagonalizability of the component matrices. 
It should be noted that a "real” version of this algorithm also exists [18]. In this latter 
algorithm, Q U Z U Q 2 , and Z 2 are orthogonal with quasi-upper triangular and R y T 
upper triangular. 
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< X < 1 


4-2. The Time-Dependent Problem: Cantilever Boundary Conditions 

A generalization of the problem which arises when modeling beams with cantilever bound- 
ary conditions is 

Cu(x,t) = !£(*,() + ^ (e/(x)|^(x,o) = /(*,(), 0 

u(o, t) = 5(0, (s/0) (1,0 = 7(0, * > 0 

^(0,0 = 3(0, k(«S|)(i.0-K0. <>o 

0 < * < 1. 


(4.16) 


du . 


u (®.°) = "^(M) = °> 

The basis for this problem is taken to be {Ct(a?)«9j(f)} where C»(®) is defined in (3.41) and 
Sj(t) = S(j, h t ) o T(t). Here the approximate solution is taken to be 

' u m.m,(®)0 = ^ 53 U *i£*{ X )Sj {t) 

i=— Af» )=-Mi 
N t 

+ J2 £?(0 { u -M m -3,j(-M a -3(x) + U-Af.-4.iC-M. -4 ( x ) 

}=- M t 

+ u Af.+i,iC^.+i(*) + Ujv.+2,iCAr.+a(*)} 

+{«(0C-M.-i(s) + ^C-M.-a(®) + 7(0Cw.+3(*) + £(<)CAr.+4(z)} 


where 


and 


7« = 


m i) 


i(t) 


FI(1)*® [JSJ(1)]» 7 ^‘ 

It should be noted that the approximate solution does satisfy the boundary conditions in 
(4.16). 

The (m* -f 4) • m t unknowns {u^} are determined by orthogonalizing the residual with 
respect to the sine functions {5 p (x)5*(t)}~ : ^.‘^^2,...,N«+2- This yields the discrete system 


(^ u m.m t ffSpS.) — 0, 


-M x -2<p<N s + 2 
—M t < q < N t 


28 


where (•, •) is defined in (4.2) with u>(x,f) — (T(t)) 

Appropriate integration by parts and the application of the one-dimensional results from 

Sections 3.2 and 3.3 yields the matrix equation 

A x UCj + C x MUA T t = G (4-17) 

where C U C X and are defined in (4.11), (4.12), and (4.13) respectively, and 

The (m* + 4) x me matrices U and ~F are defined componentwise by 

[U]ij = U ij 

and 

[F]ij tj) 

where 

7(x, t) = /(*, t) - C(a{t)B_ M .-i{ x )) - C(p{t)B_ Mm _ 2 (x)) 
-£(7(t)B w . +3 (*)) - £{&(*■) Bn.+*{ x ))- 

The (m z + 4) x (m x + 4) matrix A x is given by (3.50) or (3.52) depending on which scheme 
is used to discretize EL Finally the (m x + 4) x (m x + 4) matrix M has the form 



\ 1 

l » 

■ : 0 


M = 

a * bjj\ ! Dm 
• : 

A * A 

&H1 ! VR2 


1 ) — 

i i 0 



where the m x x m x submatrix Dm and the (m x + 4) x 1 vectors are given by 

Dm = 7}(u), 

t>L2 = D(B_m.- 4)1* 

f»Li = D(B-m.-3)1> 

&ri = D(B n .+ i)1 
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and 


b R2 = V(Bff t+2 ) 1. 

The system (4.17) can again be solved via the generalized Schur algorithm (4.15) as discussed 
in Section 4.1. 

Before implementing the method, the decay parameters 6 t / 3 , 'y and 8 must be determined 
and summation limits chosen. For the spatial weight w(x) = 1, the decay conditions are 

\EI{x)U{x,t)\ < s“ +3 (l - xf +3 t^e~ St (4.18) 

or 

\Sl{x)U{x,t)\ < *“+ 3 (l - xf +3 r + h~ st (4.19) 

depending on the manner in which El was discretized. Here SJ(x) is defined in (3.27) and 
denotes that part of the true solution which is approximated by 

«!»(*> 0 = E E u Mx)Sj ( o 

trr — Afpi J = — hft 

(see also (3.49)). With the decay parameters specified and M w chosen, the remaining stepsizes 
and summation limits are given by (4.5) - (4.9). 


5 Numerical Examples 

The four examples reported in this section were selected from a large collection of problems 
to which the Sinc-Galerkin method was applied. The results are representative of those 
obtained on other sample problems. For purposes of comparison, contrast, and performance 
evaluation, examples with known solutions were chosen. The first and last examples have 
analytic solutions, the second example has an algebraic singularity, and the third example 
contains a logarithmic singularity. As will be demonstrated by the numerical results, the 
boundary singularities have no adverse affect on the performance of the method. 
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In all examples d = *. The errors are reported on both the set of sine gridpoints 


. _ Jht 


s - {(**, x i — ^ — e< 

and the set of uniform gridpoints (l x = ~ } £ t = 

U = {(z p , J q )}p =1 , it# g 0 ; z p = £xP) 3q — ^*9 
The errors on these grids are reported as 


11^5 "(^*1 ^011 — 

-M* M\/V e 


and 


W^U ^t)ll ~ 0<p<50 u m,mi( 2: p> 

0< 9 <BO 

respectively. The dependence of both errors on the number of sine basis functions is indicated 
with the superscript M x . It is noted that if exponential convergence is realized, then 

where M x and M x denote the lower limits for the spatial sums. In the examples of this section, 
M x = 2 M X} and the exponential convergence is verified by comparing {^E^ w ^ 2 {l x , Mll)^ an< ^ 

The error and convergence results are tabulated in the form .aaa — 7 which represents 
.aaa x 10~ 7 . All problems were run with sixteen place accuracy on a Vax 8550. 


Example 5.1. 

+ 0<x<1 * >0 
u(0, <) = u(l,<) = 0, t > 0 

|(o,0-|(i.0-o, <>o 

du 

u(i,0) = —(i,0) = 0, 0 < x < 1 

The flexural rigidity is EI( x) = 1 + sin(7rx) and /(®, t) is consistent with the true solution 
u(®, t) = x(l — x) sin(47rx)t 3 e _< . The parameter is expanded via (3.25) thus yielding the 
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spatial representation (3.29). The decay condition (4.4) dictates the parameter choices ct = 
^0 = 7 = | and 8 = 1. The asymptotic error rate C^e - *''/ 3 ^/ 3 ) is maintained as indicated 
by the last column of Table 1. Notice that the choice of JV* given in (4.9) significantly 
reduces the size of the matrices involved in (4.10). The errors on both the sine grid and 
the uniform grid do not differ dramatically though in the next example the difference will 
be more noticeable. Figure 3 shows the true solution u(x, t) while Figure 4 shows both the 
true and approximate solutions (for M x = 8, 16) at the time slice t = 2. The approximate 
solution for M* = 32 is buried in the true solution on this scale. 


M t 

N x 


N t 

II 

\\&M) II 


4 

4 

4 

2 

.6207 - 0 

.8040 - 0 


8 

8 

8 

4 

.6890 - 1 

.8448 - 1 

.7345 - 0 

16 

16 

16 

6 

.2910-3 

.1105-2 

.3035 - 1 

32 

32 

32 

9 

.1906 - 4 

.2151-4 

.6587 - 4 


Table 1. Errors on the Sine Grid S and the Uniform Grid U for Example 5.1. 
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Figure 3. True Soluti 
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Example 5.2. 


+ 0 (*,<) = /(*,<), 0 <x<l 1>0 

u( 0 , t) = u(l, t) = 0 , t > 0 

^(o, 0 -go, i)-o, <>o 

«(z, o) = ^(x,0) = °> o < X < 1 

The function /(x,t) is consistent with the solution u(x,t) = [x(l — x)] 7 / a i B / a e -< which has 
algebraic singularities at x = 0, x = 1 and t — 0. The spatial discretization is taken to be 
(3.17) with the decay parameters a = /? — 7 = 2,5 = 1 dictated by (4.3). As indicated by 
Table 2 , the asymptotic error rate 0(e~' r '^“) is achieved in spite of the boundary singu- 
larities. The increased accuracy of the method for this problem as compared to Example 
5.1 is due to the larger values of a, (3 and 7 . Here the error on the sine grid is substantially 
smaller than that on the uniform grid. This emphasizes that one should use caution when 
assessing performance of a method based on only the errors at the gridpoints of the method. 
The true solution is plotted in Figure 5 while time slices (at time t = 2) of both the true 
and approximate solutions (for M x = 4, 8 ) are plotted in Figure 6 . 


M x 

N x 

M t 

N t 




4 

4 

4 

3 

.2040 - 4 

.3184-3 

" 

8 

8 

8 

4 

.9361 - 5 

.3618 - 4 

.1134-4 j 

16 

16 

16 

7 

.1372 - 6 

.1788 -5 

.5233 - 6 | 

32 

32 

32 

11 

.1742 - 9 

.1338 -7 

i 

.7441 - 8 1 


Table 2. Errors on the Sine Grid S and the Uniform Grid U for Example 5.2. 
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Example 5.3. 


d 2 u d*u 

Qi* ^ dx A ^ = 0 < aj < 1 f>0 

L 

«(°,0 = °, Q^{ M) = 2iV‘, i>0 

^(°>0 = °» * >0 

5u 

u(x, 0) = — (x,0) = 0 , 0 < x < 1 . 

The true solution u(x, t) = [(x/n(x)) 4 + x 3 ]< 3 e _i dictates the forcing function f(x,t). The 
decay condition (4.18) with EI(x ) = 1 yields the parameter choices a = fi — 1,7 = 
and 5=1 which in turn implies the asymptotic error rate (?(e -ir V ,Af "/ 2 ). As indicated by 
the results in Table 3 this rate is achieved in spite of the logarithmic singularity at x = 0. 
The convergence of the method is even accelerated which can be seen by in the last column 
of Table 3. The mesh plot in Figure 7 shows the distinctive behavior that the solution can 
exhibit when cantilever boundary conditions are in force. The “oscillation” at the right-hand 
end is tracked accurately by this method. The time slice at t = 2 shown in Figure 8 shows 
the approximate (M* = 4,8, 16) as well as the true solution. 


M x 

N x 

M t 

N t 




4 

4 

3 

2 

.4227 - 1 

.9083 - 1 


8 

8 

6 

3 

.8034 - 2 

.1999 - 1 

.3363 - 1 

16 

16 

11 

4 

.8799 - 3 

.1838 - 2 

.3953 - 2 

32 

32 

22 

7 

.3947 - 4 

.8741 -4 

.1353 - 2 


Table 3. Errors on the Sine Grid S and the Uniform Grid U for Example 5.3. 
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Example 5.4- 


d 7 u d 2 ( d 2 u \ 

^(*> 0+^1 =/(*» 0 . 0 < X < 1 <>0 

u (0,<) = 0, \EIj£j(l,t) = 0, t> 0 

1(0,0 = *-, £(*/£) (M)- 8,V«- *>0 

tt(x,0) = ^(x,0) = 0 , 0 < x < 1 . 

This example illustrates a problem where the spatial discretization (3.52) is useful. Here 
the flexural rigidity is EI(x) = 1 + sin(7rx) and the forcing function /(x, t) is consistent 
with the true solution u(x,t ) = [sin 3 (7rx) 4- x]t 2 e~‘. The parameter EI(x) is expanded via 
(3.25) thus yielding the spatial representation (3.52). The parameter choices a = (3 = 1 , 
7 =| and 8 — 1 follow from (4.19). As demonstrated by the last column of Table 4, the 
asymptotic error rate C9(e -1, V M «/ a ) j 8 achieved for larger values of M x , N X} M t , and N t . On 
this example the errors on the sine grid and those on the uniform grid are nearly the same. 
The smaller parameters a and (3 indicate why the errors here are larger than in the previous 
three examples. A mesh plot of the true solution is shown in Figure 9 while a time slice 
(< = 2) of both the true and approximate solutions ( M x = 8, 16) are plotted in Figure 10. 


M, N, M, N, \\E?-(h, ,t,)|| (l|fi" -/1 (4,<.)ll)' /5 


8 

8 

6 

3 

.6889 - 1 

.6958 - 1 


= 

16 

16 

11 

4 

.2958 - 1 

.3165-1 

.2307 - 1 

- 

32 

32 

22 

7 

.4346 - 3 

.4885 - 3 

.7572 - 2 

( 


Table 4. Errors on the Sine Grid S and the Uniform Grid U for Example 5.4. 


38 






6 Conclusions 


A fully Sinc-Galerkin method in both space and time is presented for fourth-order time- 
dependent problems with fixed and cantilever boundary conditions. The sine basis properties 
which facilitate the simple assembly of the discrete system are discussed in Section 2. In 
Section 3, the sine discretizations for the second-order temporal problem and the fourth-order 
spatial problems are presented. Alternate formulations for the variable parameter fourth- 
order problems are given which prove to be especially useful when applying the forward 
techniques of this paper to parameter recovery problems. The ODE results are then combined 
in Section 4 to form the discrete systems corresponding to the time-dependent problems of 
interest. Computational issues are discussed and a robust and efficient algorithm for solving 
the resulting matrix system is outlined. Numerical examples which highlight the method are 
given in Section 5. As demonstrated by the numerical results, the exponential convergence 
rate of the method is maintained for problems with both analytic and singular solutions as 
well as fixed and cantilever boundary conditions. 
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